function [Dtheta_Ups_tl_l_const, Dtheta_ups_tl_const] = DTheta_consts...
        (x_sub_mu_stack_s, prior, tl_gp_hypers, l_gp_hypers)
    
%    return terms required for evaluating the contribution towards the
%    variance of the bq distribution for the evidence due to marginalising
%    log input scales.

num_dims = size(x_sub_mu_stack_s, 3);

tl_input_scales = exp(tl_gp_hypers.log_input_scales);
inv_sqd_tl_input_scales_stack = ...
    reshape(tl_input_scales.^-2, 1, 1, num_dims);

L = reshape(diag(prior.covariance), 1, 1, num_dims);
Wl = reshape(exp(2*l_gp_hypers.log_input_scales), 1, 1, num_dims);
Wtl = reshape(exp(2*tl_gp_hypers.log_input_scales), 1, 1, num_dims);

denom = Wtl .* L + Wl .* L + Wl .* Wtl;
% NB: C depends on Wtl and Ct on Wl due to the inversion of a matrix
C = Wtl .* L ./ denom;
Ct = Wl .* L ./ denom;

Dtheta_Ups_tl_l_const = bsxfun(@times, inv_sqd_tl_input_scales_stack, ...
        bsxfun(@plus, ...
            L - C .* L - Ct .* L, ...
            (bsxfun(@plus, ...
                  - bsxfun(@times, C, tr(x_sub_mu_stack_s)), ... % columns
                    bsxfun(@times, (1 - Ct), x_sub_mu_stack_s))... % rows
            ).^2 ...
        )); 
                
Dtheta_ups_tl_const = bsxfun(@times, inv_sqd_tl_input_scales_stack, ...
        bsxfun(@plus, ...
            L - L .* (L + Wtl).^(-1) .* L, ...
        	bsxfun(@times, ...
                1 - L .* (L + Wtl).^(-1), ...
                x_sub_mu_stack_s).^2 ...
        ));
        